clear all
/*Insert directory here*/






/* 1) Match PUMAs */






use ctygrp70, clear
replace year=1969
append using ctygrp70
replace year=1968 if year==1970
append using ctygrp70
replace year=1967 if year==1970
append using ctygrp70
replace year=1966 if year==1970
append using ctygrp70
replace year=1965 if year==1970
append using ctygrp70
replace year=1964 if year==1970
append using ctygrp70
replace year=1963 if year==1970
append using ctygrp70
replace year=1962 if year==1970
append using ctygrp70
replace year=1961 if year==1970
append using ctygrp70
replace year=1960 if year==1970
append using ctygrp70
sort statefip geoid year
quietly by statefip geoid year: gen dup=cond(_N==1,0,_n)
drop if dup>1
drop dup
save puma2cty, replace

use ctygrp70, clear
rename year censusyear
sort statefip geoid censusyear
quietly by statefip geoid censusyear: gen dup=cond(_N==1,0,_n)
drop if dup>1
drop dup
gen furnyear=1970
save puma2ctyfurn, replace
forvalues y=1943/1969{
keep if furnyear==1970
replace furnyear=`y'
append using puma2ctyfurn
save puma2ctyfurn, replace
}

use ctygrp80, clear
replace year=1979
append using ctygrp80
replace year=1978 if year==1980
append using ctygrp80
replace year=1977 if year==1980
append using ctygrp80
replace year=1976 if year==1980
append using ctygrp80
replace year=1975 if year==1980
append using ctygrp80
replace year=1974 if year==1980
append using ctygrp80
replace year=1973 if year==1980
append using ctygrp80
replace year=1972 if year==1980
append using ctygrp80
replace year=1971 if year==1980
append using ctygrp80
sort statefip geoid year
quietly by statefip geoid year: gen dup=cond(_N==1,0,_n)
drop if dup>1
drop dup
append using puma2cty
save puma2cty, replace

use ctygrp80, clear
rename year censusyear
sort statefip geoid censusyear
quietly by statefip geoid censusyear: gen dup=cond(_N==1,0,_n)
drop if dup>1
drop dup
gen furnyear=1980
append using puma2ctyfurn
save puma2ctyfurn, replace
forvalues y=1954/1979{
keep if furnyear==1980
replace furnyear=`y'
append using puma2ctyfurn
save puma2ctyfurn, replace
}

use puma90, clear
replace year=1989
append using puma90
replace year=1988 if year==1990
append using puma90
replace year=1987 if year==1990
append using puma90
replace year=1986 if year==1990
append using puma90
replace year=1985 if year==1990
append using puma90
replace year=1984 if year==1990
append using puma90
replace year=1983 if year==1990
append using puma90
replace year=1982 if year==1990
append using puma90
replace year=1981 if year==1990
append using puma90
sort statefip geoid year
quietly by statefip geoid year: gen dup=cond(_N==1,0,_n)
drop if dup>1
drop dup
append using puma2cty
replace ctyfip=86 if ctyfip==25 & statefip==12
save puma2cty, replace

use puma90, clear
rename year censusyear
sort statefip geoid censusyear
quietly by statefip geoid censusyear: gen dup=cond(_N==1,0,_n)
drop if dup>1
drop dup
gen furnyear=1990
append using puma2ctyfurn
save puma2ctyfurn, replace
forvalues y=1964/1989{
keep if furnyear==1990
replace furnyear=`y'
append using puma2ctyfurn
save puma2ctyfurn, replace
}
replace ctyfip=86 if ctyfip==25 & statefip==12
save puma2ctyfurn, replace

use puma90-20, clear
drop if year==2020
rename puma geoid
sort statefip geoid year
quietly by statefip geoid year: gen dup=cond(_N==1,0,_n)
drop if dup>1
drop dup
append using puma2cty
save puma2cty, replace

use puma90-20, clear
drop if year<2000 | year==2020
rename puma geoid
sort statefip geoid year
quietly by statefip geoid year: gen dup=cond(_N==1,0,_n)
drop if dup>1
drop dup
rename year censusyear
gen furnyear=censusyear
append using puma2ctyfurn
save puma2ctyfurn, replace
forvalues y=1/16{
keep if furnyear==censusyear & censusyear>=2000
replace furnyear=furnyear-`y'
append using puma2ctyfurn
save puma2ctyfurn, replace
}
drop if statefip==2 | statefip==15
save puma2ctyfurn, replace

use puma2cty, clear
drop if year<1960
drop if statefip==2 | statefip==15
merge m:1 year statefip ctyfip using HDD1895-2021
keep if _merge==3
drop _merge
merge m:1 year statefip ctyfip using CDD1895-2022
keep if _merge==3
drop _merge ctyfip
save HDDgeoid60-2019, replace

use HDD1895-2021, clear
sort statefip ctyfip year
gen hddfurn=(hdd[_n-1]+hdd[_n-2]+hdd[_n-3]+hdd[_n-4]+hdd[_n-5])/5
drop hdd
rename year furnyear
merge 1:m furnyear statefip ctyfip using puma2ctyfurn
keep if _merge==3
drop _merge ctyfip
save HDDfurngeoid34-2019, replace
rename furnyear wateryear
rename hddfurn hddwater
save HDDwatergeoid34-2019, replace






/*  2) Create choice sample and run heat estimation*/






use ACS1, clear
keep if year==2000 | year==2010
drop if year==2000 & year-builtyr>10
drop if year==2010  & year-builtyr>9
gen furnyear=builtyr
rename year censusyear


/* a) Bring in outside data */


merge m:1 furnyear censusyear statefip geoid using HDDfurngeoid34-2019
keep if _merge==3
drop _merge

merge m:1 furnyear statefip using fuelpricesfurn1934-2019
keep if _merge==3
drop _merge

merge m:1 furnyear statefip censusyear geoid using incfactorfurn2019
keep if _merge==3
drop _merge
gen hhincomefurn=hhincome*incfactorfurn

merge m:1 furnyear using CPIfurn34-2019
keep if _merge==3
drop _merge

local varlist hhincomefurn pnatgasfurn poilfurn pelecfurn ppropanefurn pnatgascounterfurn
foreach a of local varlist{
replace `a'=`a'*255.7/CPI
}

replace hdd=hdd/1000
gen hdd2=hdd*hdd
gen hddfurn2=hddfurn*hddfurn

/* b) create units in structure dummies*/

gen struct2t4=0
replace struct2t4=1  if unitsstr==5 | unitsstr==6
gen struct5=0
replace struct5=1  if unitsstr>=7
gen owner=0
replace owner=1 if ownershp==1
drop ownershp unitsstr

/* e) create household size dummies*/

gen housemem2=0
replace housemem2=1 if numprec==2
gen housemem3=0
replace housemem3=1 if numprec==3
gen housemem4=0
replace housemem4=1 if numprec==4
gen housemem5=0
replace housemem5=1 if numprec==5
gen housemem6=0
replace housemem6=1 if numprec>5

/* f) HDD,rooms,members interaction terms*/

quietly sum hddfurn
gen pnatgascounterfurnxhddfurn=pnatgascounterfurn*(hddfurn-r(mean))
gen poilfurnxhddfurn=poilfurn*(hddfurn-r(mean))
gen pelecfurnxhddfurn=pelecfurn*(hddfurn-r(mean))
gen ppropanefurnxhddfurn=ppropanefurn*(hddfurn-r(mean))

quietly sum rooms
gen pnatgascounterfurnxrooms=pnatgascounterfurn*(rooms-r(mean))
gen poilfurnxrooms=poilfurn*(rooms-r(mean))
gen pelecfurnxrooms=pelecfurn*(rooms-r(mean))
gen ppropanefurnxrooms=ppropanefurn*(rooms-r(mean))

quietly sum numprec
gen pnatgascounterfurnxnumprec=pnatgascounterfurn*(numprec-r(mean))
gen poilfurnxnumprec=poilfurn*(numprec-r(mean))
gen pelecfurnxnumprec=pelecfurn*(numprec-r(mean))
gen ppropanefurnxnumprec=ppropanefurn*(numprec-r(mean))

gen systemchoice=.
replace systemchoice=1 if fuelheat==2
replace systemchoice=2 if fuelheat==5
replace systemchoice=3 if fuelheat==4
replace systemchoice=4 if fuelheat==3
save choicecensus, replace






/* 3) Run Annual Estimation Census */






clear all
/*insert directory here*/

program mylogit
  args lnf theta1 theta2 theta3 theta4 
  	 quietly replace `lnf' = ln(exp(`theta1')/(exp(`theta1')+exp(`theta2')+exp(`theta3')+exp(`theta4'))) if $ML_y1==1
     quietly replace `lnf' = ln(exp(`theta2')/(exp(`theta1')+exp(`theta2')+exp(`theta3')+exp(`theta4'))) if $ML_y1==2
	 quietly replace `lnf' = ln(exp(`theta3')/(exp(`theta1')+exp(`theta2')+exp(`theta3')+exp(`theta4'))) if $ML_y1==3
	 	 quietly replace `lnf' = ln(exp(`theta4')/(exp(`theta1')+exp(`theta2')+exp(`theta3')+exp(`theta4'))) if $ML_y1==4
end

use choicecensus, clear

constraint 1 pnatgascounterfurn=poilfurn
constraint 2 pnatgascounterfurnxhddfurn=poilfurnxhddfurn
constraint 3 pnatgascounterfurnxrooms=poilfurnxrooms
constraint 4 pnatgascounterfurnxnumprec=poilfurnxnumprec
constraint 5 ppropanefurn=poilfurn
constraint 6 ppropanefurnxhddfurn=poilfurnxhddfurn
constraint 7 ppropanefurnxrooms=poilfurnxrooms
constraint 8 ppropanefurnxnumprec=poilfurnxnumprec

local covariates hddfurn* housemem* owner rooms hhincomefurn struct*

ml model lf mylogit (gas: systemchoice=pnatgascounterfurn pnatgascounterfurnxhddfurn pnatgascounterfurnxrooms pnatgascounterfurnxnumprec `covariates')(oil: poilfurn poilfurnxhddfurn poilfurnxrooms poilfurnxnumprec, noconstant)(electricity: pelecfurn pelecfurnxhddfurn pelecfurnxrooms pelecfurnxnumprec `covariates')(propane: ppropanefurn ppropanefurnxhddfurn ppropanefurnxrooms ppropanefurnxnumprec `covariates') [pweight=hhwt], constraints (1/8)
ml check
ml search
ml maximize
clear

forvalues y=1960/2000{
use census1, clear
rename year censusyear
gen year=`y'
gen year2=year-1987
drop if censusyear>year+9 & year>1960
drop if censusyear>year+10 & year==1960
drop if censusyear<year
drop if builtyear>year

/* a) Estimate furnace and water heater replacement years */

gen unif=runiform()
gen furnyear=.
replace furnyear=builtyear if year-builtyear<=16
replace furnyear=year-0 if year-builtyear>=17 & unif>=(0/17) & unif<(1/17)
replace furnyear=year-1 if year-builtyear>=17 & unif>=(1/17) & unif<(2/17)
replace furnyear=year-2 if year-builtyear>=17 & unif>=(2/17) & unif<(3/17)
replace furnyear=year-3 if year-builtyear>=17 & unif>=(3/17) & unif<(4/17)
replace furnyear=year-4 if year-builtyear>=17 & unif>=(4/17) & unif<(5/17)
replace furnyear=year-5 if year-builtyear>=17 & unif>=(5/17) & unif<(6/17)
replace furnyear=year-6 if year-builtyear>=17 & unif>=(6/17) & unif<(7/17)
replace furnyear=year-7 if year-builtyear>=17 & unif>=(7/17) & unif<(8/17)
replace furnyear=year-8 if year-builtyear>=17 & unif>=(8/17) & unif<(9/17)
replace furnyear=year-9 if year-builtyear>=17 & unif>=(9/17) & unif<(10/17)
replace furnyear=year-10 if year-builtyear>=17 & unif>=(10/17) & unif<(11/17)
replace furnyear=year-11 if year-builtyear>=17 & unif>=(11/17) & unif<(12/17)
replace furnyear=year-12 if year-builtyear>=17 & unif>=(12/17) & unif<(13/17)
replace furnyear=year-13 if year-builtyear>=17 & unif>=(13/17) & unif<(14/17)
replace furnyear=year-14 if year-builtyear>=17 & unif>=(14/17) & unif<(15/17)
replace furnyear=year-15 if year-builtyear>=17 & unif>=(15/17) & unif<(16/17)
replace furnyear=year-16 if year-builtyear>=17 & unif>=(16/17) & unif<=(17/17)

gen wateryear=.
replace wateryear=builtyear if year-builtyear<=9
replace wateryear=year-0 if year-builtyear>=10 & unif>=(0/10) & unif<(1/10)
replace wateryear=year-1 if year-builtyear>=10 & unif>=(1/10) & unif<(2/10)
replace wateryear=year-2 if year-builtyear>=10 & unif>=(2/10) & unif<(3/10)
replace wateryear=year-3 if year-builtyear>=10 & unif>=(3/10) & unif<(4/10)
replace wateryear=year-4 if year-builtyear>=10 & unif>=(4/10) & unif<(5/10)
replace wateryear=year-5 if year-builtyear>=10 & unif>=(5/10) & unif<(6/10)
replace wateryear=year-6 if year-builtyear>=10 & unif>=(6/10) & unif<(7/10)
replace wateryear=year-7 if year-builtyear>=10 & unif>=(7/10) & unif<(8/10)
replace wateryear=year-8 if year-builtyear>=10 & unif>=(8/10) & unif<(9/10)
replace wateryear=year-9 if year-builtyear>=10 & unif>=(9/10) & unif<=(10/10)
drop unif

/* b) Bring in outside data */

/*current year*/
merge m:1 year statefip geoid using HDDgeoid60-2019
keep if _merge==3
drop _merge

merge m:1 year statefip using incfactor2000
keep if _merge==3
drop _merge
replace hhincome=hhincome*incfactor

merge m:1 year statefip using fuelprices1960-19
keep if _merge==3
drop _merge

merge m:1 year using CPI34-2019
keep if _merge==3
drop _merge

/*furnace year*/

merge m:1 furnyear censusyear statefip geoid using HDDfurngeoid34-2019
keep if _merge==3
drop _merge

merge m:1 furnyear statefip censusyear using incfactorfurn2000
keep if _merge==3
drop _merge
gen hhincomefurn=hhincome*incfactorfurn

merge m:1 furnyear statefip using fuelpricesfurn1934-2019
keep if _merge==3
drop _merge

merge m:1 furnyear using CPIfurn34-2019
keep if _merge==3
drop _merge

/*water year*/

merge m:1 wateryear censusyear statefip geoid using HDDwatergeoid34-2019
keep if _merge==3
drop _merge

merge m:1 wateryear statefip censusyear using incfactorwater2000
keep if _merge==3
drop _merge
gen hhincomewater=hhincome*incfactorwater

merge m:1 wateryear statefip using fuelpriceswater1934-2019
keep if _merge==3
drop _merge

merge m:1 wateryear using CPIwater34-2019
keep if _merge==3
drop _merge

local varlist hhincome pnatgas poil pelec ppropane pnatgascounter
foreach b of local varlist{
replace `b'=`b'*255.7/CPI
}

local varlist2 hhincomefurn pnatgasfurn poilfurn pelecfurn ppropanefurn pnatgascounterfurn
foreach c of local varlist2{
replace `c'=`c'*255.7/CPIfurn
}

local varlist3 hhincomewater pnatgaswater poilwater pelecwater ppropanewater pnatgascounterwater
foreach d of local varlist3{
replace `d'=`d'*255.7/CPIwater
}

drop incfactor incfactorfurn incfactorwater
replace hdd=hdd/1000
replace hddfurn=hddfurn/1000
replace hddwater=hddwater/1000
gen hdd2=hdd*hdd
gen hddfurn2=hddfurn*hddfurn
gen hddwater2=hddwater*hddwater
replace cdd=cdd/1000
gen cdd2=cdd*cdd

/* c) create units in structure dummies*/

gen struct2t4=0
replace struct2t4=1  if unitsstr==5 | unitsstr==6
gen struct5=0
replace struct5=1  if unitsstr>=7
gen owner=0
replace owner=1 if ownershp==1
drop ownershp unitsstr

/* d) create age of structure dummies*/

gen ageb1950s=0
gen ageb1960s=0
gen ageb1970s=0
gen ageb1980s=0
gen ageb1990s=0
gen ageb2000s=0
drop if builtyr==1 & censusyear==1960 & year<1959
drop if builtyr==2 & censusyear==1960 & year<1955
replace ageb1950s=1 & censusyear==1960 & builtyr<=3

drop if builtyr==1 & censusyear==1970 & year<1969
drop if builtyr==2 & censusyear==1970 & year<1965
replace ageb1960s=1 if censusyear==1970 & builtyr<=3
replace ageb1950s=1 if censusyear==1970 & builtyr==4

drop if builtyr==1 & censusyear==1980 & year<1979
drop if builtyr==2 & censusyear==1980 & year<1975
replace ageb1970s=1 if censusyear==1980 & builtyr<=3
replace ageb1960s=1 if censusyear==1980 & builtyr==4
replace ageb1950s=1 if censusyear==1980 & builtyr==5

drop if builtyr==1 & censusyear==1990 & year<1989
drop if builtyr==2 & censusyear==1990 & year<1985
replace ageb1980s=1 if censusyear==1990 & builtyr<=3
replace ageb1970s=1 if censusyear==1990 & builtyr==4
replace ageb1960s=1 if censusyear==1990 & builtyr==5
replace ageb1950s=1 if censusyear==1990 & builtyr==6

drop if builtyr==1 & censusyear==2000 & year<1999
drop if builtyr==2 & censusyear==2000 & year<1995
replace ageb1990s=1 if censusyear==2000 & builtyr<=3
replace ageb1980s=1 if censusyear==2000 & builtyr==4
replace ageb1970s=1 if censusyear==2000 & builtyr==5
replace ageb1960s=1 if censusyear==2000 & builtyr==6
replace ageb1950s=1 if censusyear==2000 & builtyr==7

/* e) create household size dummies*/

gen housemem2=0
replace housemem2=1 if numprec==2
gen housemem3=0
replace housemem3=1 if numprec==3
gen housemem4=0
replace housemem4=1 if numprec==4
gen housemem5=0
replace housemem5=1 if numprec==5
gen housemem6=0
replace housemem6=1 if numprec>5

/* f) HDD,rooms,members interaction terms*/

local varlist furn water
foreach i of local varlist{
quietly sum hdd`i'
gen pnatgascounter`i'xhdd`i'=pnatgascounter`i'*(hdd`i'-r(mean))
gen poil`i'xhdd`i'=poil`i'*(hdd`i'-r(mean))
gen pelec`i'xhdd`i'=pelec`i'*(hdd`i'-r(mean))
gen ppropane`i'xhdd`i'=ppropane`i'*(hdd`i'-r(mean))

quietly sum rooms
gen pnatgascounter`i'xrooms=pnatgascounter`i'*(rooms-r(mean))
gen poil`i'xrooms=poil`i'*(rooms-r(mean))
gen pelec`i'xrooms=pelec`i'*(rooms-r(mean))
gen ppropane`i'xrooms=ppropane`i'*(rooms-r(mean))

quietly sum numprec
gen pnatgascounter`i'xnumprec=pnatgascounter`i'*(numprec-r(mean))
gen poil`i'xnumprec=poil`i'*(numprec-r(mean))
gen pelec`i'xnumprec=pelec`i'*(numprec-r(mean))
gen ppropane`i'xnumprec=ppropane`i'*(numprec-r(mean))
}

/* g) Create predicted values of space heating fuel choice */

gen predg=[gas]_b[pnatgascounterfurn]*pnatgascounterfurn+[gas]_b[pnatgascounterfurnxhddfurn]*pnatgascounterfurnxhddfurn+[gas]_b[pnatgascounterfurnxrooms]*pnatgascounterfurnxrooms+[gas]_b[pnatgascounterfurnxnumprec]*pnatgascounterfurnxnumprec+[gas]_b[hddfurn]*hddfurn+[gas]_b[hddfurn2]*hddfurn2+[gas]_b[housemem2]*housemem2+[gas]_b[housemem3]*housemem3+[gas]_b[housemem4]*housemem4+[gas]_b[housemem5]*housemem5+[gas]_b[housemem6]*housemem6+[gas]_b[owner]*owner+[gas]_b[hhincome]*hhincomefurn+[gas]_b[rooms]*rooms+[gas]_b[struct2t4]*struct2t4+[gas]_b[struct5]*struct5+[gas]_b[_cons]

gen predo=[oil]_b[poilfurn]*poilfurn+[oil]_b[poilfurnxhddfurn]*poilfurnxhddfurn+[oil]_b[poilfurnxrooms]*poilfurnxrooms+[oil]_b[poilfurnxnumprec]*poilfurnxnumprec

gen prede=[electricity]_b[pelecfurn]*pelecfurn+[electricity]_b[pelecfurnxhddfurn]*pelecfurnxhddfurn+[electricity]_b[pelecfurnxrooms]*pelecfurnxrooms+[electricity]_b[pelecfurnxnumprec]*pelecfurnxnumprec+[electricity]_b[hddfurn]*hddfurn+[electricity]_b[hddfurn2]*hddfurn2+[electricity]_b[housemem2]*housemem2+[electricity]_b[housemem3]*housemem3+[electricity]_b[housemem4]*housemem4+[electricity]_b[housemem5]*housemem5+[electricity]_b[housemem6]*housemem6+[electricity]_b[owner]*owner+[electricity]_b[hhincome]*hhincomefurn+[electricity]_b[rooms]*rooms+[electricity]_b[struct2t4]*struct2t4+[electricity]_b[struct5]*struct5+[electricity]_b[_cons]

gen predp=[propane]_b[ppropanefurn]*ppropanefurn+[propane]_b[ppropanefurnxhddfurn]*ppropanefurnxhddfurn+[propane]_b[ppropanefurnxrooms]*ppropanefurnxrooms+[propane]_b[ppropanefurnxnumprec]*ppropanefurnxnumprec+[propane]_b[hddfurn]*hddfurn+[propane]_b[hddfurn2]*hddfurn2+[propane]_b[housemem2]*housemem2+[propane]_b[housemem3]*housemem3+[propane]_b[housemem4]*housemem4+[propane]_b[housemem5]*housemem5+[propane]_b[housemem6]*housemem6+[propane]_b[owner]*owner+[propane]_b[hhincome]*hhincomefurn+[propane]_b[rooms]*rooms+[propane]_b[struct2t4]*struct2t4+[propane]_b[struct5]*struct5+[propane]_b[_cons]

gen choiceg=exp(predg)/(exp(predg)+exp(predo)+exp(prede)+exp(predp))
gen choiceo=exp(predo)/(exp(predg)+exp(predo)+exp(prede)+exp(predp))
gen choicee=exp(prede)/(exp(predg)+exp(predo)+exp(prede)+exp(predp))
gen choicep=exp(predp)/(exp(predg)+exp(predo)+exp(prede)+exp(predp))

save dataset`y', replace
}






/* 4) Run RECS water heater choice model for Census */






clear all
program mylogit
  args lnf theta1 theta2
  	 quietly replace `lnf' = ln(exp(`theta1')/(exp(`theta1')+exp(`theta2'))) if $ML_y1==1
     quietly replace `lnf' = ln(exp(`theta2')/(exp(`theta1')+exp(`theta2'))) if $ML_y1==2
end


/* a) Water heater with gas furnace */


use RECSadj, clear
keep if fuelheat==1
keep if fuelh2o==1 | fuelh2o==5

gen systemchoice=.
replace systemchoice=1 if fuelh2o==5
replace systemchoice=2 if fuelh2o==1

local covariates hdd* owner rooms hhincome struct* housemem*

ml model lf mylogit (elec: systemchoice=elecprice elecpricexhdd elecpricexrooms elecpricexhhmem  `covariates')(gas: gasprice gaspricexhdd gaspricexrooms gaspricexhhmem, noconstant) [pweight=nweight]
ml check
ml search
ml maximize
clear

forvalues y=1960/2000{
use dataset`y', clear

/* electric with gas furnace */

gen predwaterge=[elec]_b[elecprice]*pelecwater+[elec]_b[elecpricexhdd]*pelecwaterxhddwater+[elec]_b[elecpricexrooms]*pelecwaterxrooms+[elec]_b[elecpricexhhmem]*pelecwaterxnumprec+[elec]_b[hdd]*hddwater+[elec]_b[hdd2]*hddwater2+[elec]_b[housemem2]*housemem2+[elec]_b[housemem3]*housemem3+[elec]_b[housemem4]*housemem4+[elec]_b[housemem5]*housemem5+[elec]_b[housemem6]*housemem6+[elec]_b[owner]*owner+[elec]_b[hhincome]*hhincome+[elec]_b[rooms]*rooms+[elec]_b[struct2t4]*struct2t4+[elec]_b[struct5]*struct5+[elec]_b[_cons]

/* gas with gas furnace */
gen predwatergg=[gas]_b[gasprice]*pnatgascounterwater+[gas]_b[gaspricexhdd]*pnatgascounterwaterxhddwater+[gas]_b[gaspricexrooms]*pnatgascounterwaterxrooms+[gas]_b[gaspricexhhmem]*pnatgascounterwaterxnumprec

gen choicewaterge=exp(predwaterge)/(exp(predwaterge)+exp(predwatergg))
gen choicewatergg=exp(predwatergg)/(exp(predwaterge)+exp(predwatergg))
drop predwaterge predwatergg
save dataset`y', replace
}


/* b) Water heater with electric furnace */


use RECSadj, clear
keep if fuelheat==5
keep if fuelh2o==1 | fuelh2o==5

gen systemchoice=.
replace systemchoice=1 if fuelh2o==1
replace systemchoice=2 if fuelh2o==5

local covariates hdd* owner rooms hhincome struct* housemem*

ml model lf mylogit (gas: systemchoice=gasprice gaspricexhdd gaspricexrooms gaspricexhhmem  `covariates')(elec: elecprice elecpricexhdd elecpricexrooms elecpricexhhmem, noconstant) [pweight=nweight]
ml check
ml search
ml maximize
clear

forvalues y=1960/2000{
use dataset`y', clear

/* gas with electric furnace */
gen predwatereg=[gas]_b[gasprice]*pnatgascounterwater+[gas]_b[gaspricexhdd]*pnatgascounterwaterxhddwater+[gas]_b[gaspricexrooms]*pnatgascounterwaterxrooms+[gas]_b[gaspricexhhmem]*pnatgascounterwaterxnumprec+[gas]_b[hdd]*hddwater+[gas]_b[hdd2]*hddwater2+[gas]_b[housemem2]*housemem2+[gas]_b[housemem3]*housemem3+[gas]_b[housemem4]*housemem4+[gas]_b[housemem5]*housemem5+[gas]_b[housemem6]*housemem6+[gas]_b[owner]*owner+[gas]_b[hhincome]*hhincome+[gas]_b[rooms]*rooms+[gas]_b[struct2t4]*struct2t4+[gas]_b[struct5]*struct5+[gas]_b[_cons]

/*electric with electric furnace*/
gen predwateree=[elec]_b[elecprice]*pelecwater+[elec]_b[elecpricexhdd]*pelecwaterxhddwater+[elec]_b[elecpricexrooms]*pelecwaterxrooms+[elec]_b[elecpricexhhmem]*pelecwaterxnumprec

gen choicewatereg=exp(predwatereg)/(exp(predwatereg)+exp(predwateree))
gen choicewateree=exp(predwateree)/(exp(predwatereg)+exp(predwateree))
drop predwateree predwatereg
save dataset`y', replace
}


/* c) Water heater with oil furnace*/


use RECSadj, clear
keep if fuelheat==3
keep if fuelh2o==3 | fuelh2o==5

gen systemchoice=.
replace systemchoice=1 if fuelh2o==5
replace systemchoice=2 if fuelh2o==3

local covariates hdd* owner rooms hhincome struct* housemem*

ml model lf mylogit (elec: systemchoice=elecprice elecpricexhdd elecpricexrooms elecpricexhhmem  `covariates')(oil: oilprice oilpricexhdd oilpricexrooms oilpricexhhmem, noconstant) [pweight=nweight]
ml check
ml search
ml maximize
clear

forvalues y=1960/2000{
use dataset`y', clear

/*electric with oil furnace*/
gen predwateroe=[elec]_b[elecprice]*pelecwater+[elec]_b[elecpricexhdd]*pelecwaterxhddwater+[elec]_b[elecpricexrooms]*pelecwaterxrooms+[elec]_b[elecpricexhhmem]*pelecwaterxnumprec+[elec]_b[hdd]*hddwater+[elec]_b[hdd2]*hddwater2+[elec]_b[housemem2]*housemem2+[elec]_b[housemem3]*housemem3+[elec]_b[housemem4]*housemem4+[elec]_b[housemem5]*housemem5+[elec]_b[housemem6]*housemem6+[elec]_b[owner]*owner+[elec]_b[hhincome]*hhincome+[elec]_b[rooms]*rooms+[elec]_b[struct2t4]*struct2t4+[elec]_b[struct5]*struct5+[elec]_b[_cons]

/*oil with oil furnace*/
gen predwateroo=[oil]_b[oilprice]*poilwater+[oil]_b[oilpricexhdd]*poilwaterxhddwater+[oil]_b[oilpricexrooms]*poilwaterxrooms+[oil]_b[oilpricexhhmem]*poilwaterxnumprec

gen choicewateroe=exp(predwateroe)/(exp(predwateroe)+exp(predwateroo))
gen choicewateroo=exp(predwateroo)/(exp(predwateroe)+exp(predwateroo))
drop predwateroo predwateroe
save dataset`y', replace
}


/* d) Propane furnace water heater fuel */

use RECSadj, clear
keep if fuelheat==2
keep if fuelh2o==2 | fuelh2o==5

gen systemchoice=.
replace systemchoice=1 if fuelh2o==5
replace systemchoice=2 if fuelh2o==2

local covariates hdd* owner rooms hhincome struct* housemem*

ml model lf mylogit (elec: systemchoice=elecprice elecpricexhdd elecpricexrooms elecpricexhhmem  `covariates')(prop: propprice proppricexhdd proppricexrooms proppricexhhmem, noconstant) [pweight=nweight]
ml check
ml search
ml maximize
clear

forvalues y=1960/2000{
use dataset`y', clear

/*electric with prop furnace*/
gen predwaterpe=[elec]_b[elecprice]*pelecwater+[elec]_b[elecpricexhdd]*pelecwaterxhddwater+[elec]_b[elecpricexrooms]*pelecwaterxrooms+[elec]_b[elecpricexhhmem]*pelecwaterxnumprec+[elec]_b[hdd]*hddwater+[elec]_b[hdd2]*hddwater2+[elec]_b[housemem2]*housemem2+[elec]_b[housemem3]*housemem3+[elec]_b[housemem4]*housemem4+[elec]_b[housemem5]*housemem5+[elec]_b[housemem6]*housemem6+[elec]_b[owner]*owner+[elec]_b[hhincome]*hhincome+[elec]_b[rooms]*rooms+[elec]_b[struct2t4]*struct2t4+[elec]_b[struct5]*struct5+[elec]_b[_cons]

/*prop with prop furnace*/
gen predwaterpp=[prop]_b[propprice]*ppropanewater+[prop]_b[proppricexhdd]*ppropanewaterxhddwater+[prop]_b[proppricexrooms]*ppropanewaterxrooms+[prop]_b[proppricexhhmem]*ppropanewaterxnumprec

gen choicewaterpe=exp(predwaterpe)/(exp(predwaterpe)+exp(predwaterpp))
gen choicewaterpp=exp(predwaterpp)/(exp(predwaterpe)+exp(predwaterpp))
drop predwaterpp predwaterpe

gen choicewaterg=choiceg*choicewatergg+choicee*choicewatereg
gen choicewatero=choiceo*choicewateroo
gen choicewaterp=choicep*choicewaterpp
gen choicewatere=choicee*choicewateree+choiceg*choicewaterge+choiceo*choicewateroe+choicep*choicewaterpe
save dataset`y', replace
}






/* 5) Generate Demand Census*/






/* a) Electric Space Heating Demand*/

forvalues y=1960/1980{

use RECSadj, clear
gen id1=0
replace id1=1 if year<=1993 & year>=1987
gen id2=0
replace id2=1 if year<=1997 & year>=1987
gen id3=0
replace id3=1 if year<=2001 & year>=1987

drop if id1!=1 & `y'<1980
drop if id2!=1 & `y'>=1980 & `y'<1990
drop if id3!=1 & `y'>=1990
drop id1 id2 id3

keep if fuelheat==5
reg heatcons heatprice hdd hdd2 hhincome rooms build* struct* owner [pweight=nweight], robust

use dataset`y', clear
drop predelecspacecons prednatgasspacecons predoilspacecons predpropanespacecons predelecwatercons predpropanewatercons predoilwatercons prednatgaswatercons predcomeleccons predac

gen predelecspacecons=_b[heatprice]*pelec+_b[hdd]*hdd+_b[hdd2]*hdd2+_b[owner]*owner+_b[hhincome]*hhincome+_b[rooms]*rooms+_b[struct2t4]*struct2t4+_b[struct5]*struct5+_b[build1950s]*ageb1950+_b[build1960s]*ageb1960+_b[build1970s]*ageb1970+_b[build1980s]*ageb1980+_b[build1990s]*ageb1990+_b[_cons]
replace predelecspacecons=0 if predelecspacecons<0

save dataset`y', replace
}

forvalues y=1981/2000{

use RECSadj, clear
gen id1=0
replace id1=1 if year<=1993 & year>=1987
gen id2=0
replace id2=1 if year<=1997 & year>=1987
gen id3=0
replace id3=1 if year<=2001 & year>=1987

drop if id1!=1 & `y'<1980
drop if id2!=1 & `y'>=1981 & `y'<1990
drop if id3!=1 & `y'>=1990
drop id1 id2 id3

keep if fuelheat==5
reg heatcons heatprice hdd hdd2 hhincome rooms build* struct* owner [pweight=nweight], robust

use dataset`y', clear
drop predelecspacecons prednatgasspacecons predoilspacecons predpropanespacecons predelecwatercons predpropanewatercons predoilwatercons prednatgaswatercons predcomeleccons

gen predelecspacecons=_b[heatprice]*pelec+_b[hdd]*hdd+_b[hdd2]*hdd2+_b[owner]*owner+_b[hhincome]*hhincome+_b[rooms]*rooms+_b[struct2t4]*struct2t4+_b[struct5]*struct5+_b[build1950s]*ageb1950+_b[build1960s]*ageb1960+_b[build1970s]*ageb1970+_b[build1980s]*ageb1980+_b[build1990s]*ageb1990+_b[_cons]
replace predelecspacecons=0 if predelecspacecons<0

save dataset`y', replace
}

/* b) NatGas & Oil Space Heating Demand*/

forvalues y=1960/1989{
use RECSadj, clear
keep if year>=1987 & year<=1993
keep if fuelheat==1 | fuelheat==3
reg heatcons heatprice hdd hdd2 hhincome rooms build* struct* owner year2 [pweight=nweight], robust

use dataset`y', clear

rename pnatgascounter pnatgas
local varlist natgas oil
foreach j of local varlist{
gen pred`j'spacecons=_b[heatprice]*p`j'+_b[hdd]*hdd+_b[hdd2]*hdd2+_b[owner]*owner+_b[hhincome]*hhincome+_b[rooms]*rooms+_b[struct2t4]*struct2t4+_b[struct5]*struct5+_b[build1950s]*ageb1950+_b[build1960s]*ageb1960+_b[year2]*year2+_b[_cons]
replace pred`j'spacecons=pred`j'spacecons+_b[build1970s]*ageb1970 if year>=1970
replace pred`j'spacecons=pred`j'spacecons+_b[build1980s]*ageb1980 if year>=1980
replace pred`j'spacecons=0 if pred`j'spacecons<0
}
rename pnatgas pnatgascounter
save dataset`y', replace
}


forvalues y=1990/2000{
use RECSadj, clear
gen id1=0
replace id1=1 if year==1993
gen id2=0
replace id2=1 if year>=1997 & year<=2001
gen id3=0
replace id3=1 if year>=1997 & year<=2005

drop if id1!=1 & `y'>=1990 & `y'<1994
drop if id2!=1 & `y'>=1994 & `y'<1997
drop if id3!=1 & `y'>=1997 & `y'<=2000
drop id1 id2 id3
keep if fuelheat==1 | fuelheat==3
reg heatcons heatprice hdd hdd2 hhincome rooms build* struct* owner [pweight=nweight], robust

use dataset`y', clear

rename pnatgascounter pnatgas
local varlist natgas oil
foreach j of local varlist{
gen pred`j'spacecons=_b[heatprice]*p`j'+_b[hdd]*hdd+_b[hdd2]*hdd2+_b[owner]*owner+_b[hhincome]*hhincome+_b[rooms]*rooms+_b[struct2t4]*struct2t4+_b[struct5]*struct5+_b[build1950s]*ageb1950+_b[build1960s]*ageb1960+_b[_cons]
replace pred`j'spacecons=pred`j'spacecons+_b[build1970s]*ageb1970 if year>=1970
replace pred`j'spacecons=pred`j'spacecons+_b[build1980s]*ageb1980 if year>=1980
replace pred`j'spacecons=pred`j'spacecons+_b[build1990s]*ageb1990 if year>=1990
replace pred`j'spacecons=0 if pred`j'spacecons<0
}
rename pnatgas pnatgascounter
save dataset`y', replace
}

/* c) Propane Space Heating Demand*/

forvalues y=1960/2000{
use RECSadj, clear
gen id1=0
replace id1=1 if year==1987
gen id2=0
replace id2=1 if year==1993
gen id3=0
replace id3=1 if year==1997
gen id4=0
replace id4=1 if year>=1997 & year<=2005

drop if id1!=1 & `y'<1966
drop if id2!=1 & `y'>=1966 & `y'<1979
drop if id3!=1 & `y'>=1979 & `y'<1990
drop if id4!=1 & `y'>=1990 & `y'<=2000
drop id1 id2 id3 id4
keep if fuelheat==2
reg heatcons heatprice hdd hdd2 hhincome rooms build* struct* owner [pweight=nweight], robust

use dataset`y', clear

gen predpropanespacecons=_b[heatprice]*ppropane+_b[hdd]*hdd+_b[hdd2]*hdd2+_b[owner]*owner+_b[hhincome]*hhincome+_b[rooms]*rooms+_b[struct2t4]*struct2t4+_b[struct5]*struct5+_b[build1950s]*ageb1950+_b[build1960s]*ageb1960+_b[_cons]
replace predpropanespacecons=predpropanespacecons+_b[build1970s]*ageb1970 if year>=1970
replace predpropanespacecons=predpropanespacecons+_b[build1980s]*ageb1980 if year>=1980
replace predpropanespacecons=predpropanespacecons+_b[build1990s]*ageb1990 if year>=1990
replace predpropanespacecons=0 if predpropanespacecons<0

save dataset`y', replace
}

/* d) Generate Water Heating Demand */


/*Electric Water Heating*/

forvalues y=1960/2000{
use RECSadj, clear
gen id1=0
replace id1=1 if year<=1993 & year>=1987
gen id2=0
replace id2=1 if year<=1997 & year>=1987
gen id3=0
replace id3=1 if year<=2001 & year>=1987

drop if id1!=1 & `y'<1980
drop if id2!=1 & `y'>=1980 & `y'<1990
drop if id3!=1 & `y'>=1990
drop id1 id2 id3

keep if fuelh2o==5
drop if heatprice<1 | heatprice>300
drop if watercons<0.5 | watercons>80

reg watercons elecprice hdd hdd2 hhincome rooms struct* housemem* owner [pweight=nweight], robust

use dataset`y', clear

gen predelecwatercons=_b[elecprice]*pelec+_b[hdd]*hdd+_b[hdd2]*hdd2+_b[owner]*owner+_b[hhincome]*hhincome+_b[rooms]*rooms+_b[struct2t4]*struct2t4+_b[struct5]*struct5+_b[_cons]
replace predelecwatercons=0 if predelecwatercons<0

save dataset`y', replace
}

/*Fossil Fuel Water Heating*/

forvalues y=1960/2000{

use RECSadj, clear
gen id1=0
replace id1=1 if year<=1993 & year>=1987
gen id2=0
replace id2=1 if year<=1997 & year>=1987
gen id3=0
replace id3=1 if year<=2001 & year>=1987

drop if id1!=1 & `y'<1980
drop if id2!=1 & `y'>=1980 & `y'<1990
drop if id3!=1 & `y'>=1990
drop id1 id2 id3

drop if fuelh2o==5
drop if heatprice<1 | heatprice>300
drop if watercons<0.5 | watercons>80

reg watercons waterheatprice hdd hdd2 hhincome rooms struct* housemem* owner [pweight=nweight], robust

use dataset`y', clear

local varlist2 propane oil
foreach j of local varlist2{
gen pred`j'watercons=_b[waterheatprice]*p`j'+_b[hdd]*hdd+_b[hdd2]*hdd2+_b[owner]*owner+_b[hhincome]*hhincome+_b[rooms]*rooms+_b[struct2t4]*struct2t4+_b[struct5]*struct5+_b[_cons]
replace pred`j'watercons=0 if pred`j'watercons<0
}

gen prednatgaswatercons=_b[waterheatprice]*pnatgascounter+_b[hdd]*hdd+_b[hdd2]*hdd2+_b[owner]*owner+_b[hhincome]*hhincome+_b[rooms]*rooms+_b[struct2t4]*struct2t4+_b[struct5]*struct5+_b[_cons]
replace prednatgaswatercons=0 if prednatgaswatercons<0

save dataset`y', replace
}


/* e) Generate Non-heating electricity demand */


use dataset1970, clear
append using dataset1980

reg aircentral pelec cdd cdd2 hhincome rooms ageb1950s ageb1960s ageb1970s struct* owner [pweight=hhwt], robust

forvalues y=1960/1980{
use dataset`y', clear

gen predac=_b[pelec]*pelec+_b[cdd]*cdd+_b[cdd2]*cdd2+_b[owner]*owner+_b[hhincome]*hhincome+_b[rooms]*rooms+_b[struct2t4]*struct2t4+_b[struct5]*struct5+_b[ageb1950s]*ageb1950+_b[ageb1960s]*ageb1960+_b[ageb1970s]*ageb1970+_b[_cons]
replace predac=0 if predac<0
replace predac=1 if predac>1
save dataset`y', replace
}

forvalues y=1960/1980{
use RECSadj, clear

gen id1=0
replace id1=1 if year<1987
gen id2=0
replace id2=1 if year<=1987
gen id3=0
replace id3=1 if year<=1990
gen id4=0
replace id4=1 if year<=1997
gen id5=0
replace id5=1 if year<=2001


drop if id1!=1 & `y'<1961
drop if id2!=1 & `y'>=1961 & `y'<1969
drop if id3!=1 & `y'>=1969 & `y'<1974
drop if id4!=1 & `y'>=1974 & `y'<1977
drop if id5!=1 & `y'>=1977
drop id1 id2 id3 id4 id5

reg comeleccons elecprice cdd cdd2 hhincome rooms build* struct* owner aircentral [pweight=nweight], robust


use dataset`y', clear

gen predcomeleccons=_b[elecprice]*pelec+_b[cdd]*cdd+_b[cdd2]*cdd2+_b[owner]*owner+_b[hhincome]*hhincome+_b[rooms]*rooms+_b[struct2t4]*struct2t4+_b[struct5]*struct5+_b[build1950s]*ageb1950+_b[build1960s]*ageb1960+_b[build1970s]*ageb1970+_b[build1980s]*ageb1980+_b[build1990s]*ageb1990+_b[aircentral]*predac+_b[_cons]
replace predcomeleccons=0 if predcomeleccons<0

save dataset`y', replace
}

forvalues y=1981/2000{

use RECSadj, clear

gen id1=0
replace id1=1 if year<=1990 & year>=1987
gen id2=0
replace id2=1 if year<=1997 & year>=1990
gen id3=0
replace id3=1 if year<=2001 & year>=1993
gen id4=0
replace id4=1 if year<=2001 & year>=1997

drop if id1!=1 & `y'>1980 & `y'<1987
drop if id2!=1 & `y'>=1987 & `y'<1993
drop if id3!=1 & `y'>=1993 & `y'<1998
drop if id4!=1 & `y'>=1998
drop id1 id2 id3 id4

reg comeleccons elecprice cdd cdd2 hhincome rooms build* struct* owner aircentral [pweight=nweight], robust


use dataset`y', clear

gen predcomeleccons=_b[elecprice]*pelec+_b[cdd]*cdd+_b[cdd2]*cdd2+_b[owner]*owner+_b[hhincome]*hhincome+_b[rooms]*rooms+_b[struct2t4]*struct2t4+_b[struct5]*struct5+_b[build1950s]*ageb1950+_b[build1960s]*ageb1960+_b[build1970s]*ageb1970+_b[build1980s]*ageb1980+_b[build1990s]*ageb1990+_b[_cons]
replace predcomeleccons=0 if predcomeleccons<0

save dataset`y', replace

}






/* 6) Create county averages */






/*Convert geoids to counties*/

use puma90-20, clear
sort year statefip ctyfip
quietly by year statefip ctyfip: gen dup=cond(_N==1,0,_n)
drop if dup>1
drop dup
rename puma geoid
save geoid2cty, replace

use puma90, clear
replace year=1989
append using puma90
replace year=1988 if year==1990
append using puma90
replace year=1987 if year==1990
append using puma90
replace year=1986 if year==1990
append using puma90
replace year=1985 if year==1990
append using puma90
replace year=1984 if year==1990
append using puma90
replace year=1983 if year==1990
append using puma90
replace year=1982 if year==1990
append using puma90
replace year=1981 if year==1990
append using puma90
sort year statefip ctyfip
quietly by year statefip ctyfip: gen dup=cond(_N==1,0,_n)
drop if dup>1
drop dup
append using geoid2cty
save geoid2cty, replace

use ctygrp80, clear
replace year=1979
append using ctygrp80
replace year=1978 if year==1980
append using ctygrp80
replace year=1977 if year==1980
append using ctygrp80
replace year=1976 if year==1980
append using ctygrp80
replace year=1975 if year==1980
append using ctygrp80
replace year=1974 if year==1980
append using ctygrp80
replace year=1973 if year==1980
append using ctygrp80
replace year=1972 if year==1980
append using ctygrp80
replace year=1971 if year==1980
append using ctygrp80
sort year statefip ctyfip
quietly by year statefip ctyfip: gen dup=cond(_N==1,0,_n)
drop if dup>1
drop dup
append using geoid2cty
save geoid2cty, replace

use ctygrp70, clear
replace year=1969
append using ctygrp70
replace year=1968 if year==1970
append using ctygrp70
replace year=1967 if year==1970
append using ctygrp70
replace year=1966 if year==1970
append using ctygrp70
replace year=1965 if year==1970
append using ctygrp70
replace year=1964 if year==1970
append using ctygrp70
replace year=1963 if year==1970
append using ctygrp70
replace year=1962 if year==1970
append using ctygrp70
replace year=1961 if year==1970
append using ctygrp70
replace year=1960 if year==1970
append using ctygrp70
sort year statefip ctyfip
quietly by year statefip ctyfip: gen dup=cond(_N==1,0,_n)
drop if dup>1
drop dup
append using geoid2cty
drop if statefip==2 | statefip==15
save geoid2cty, replace

/*Generate average by geoid*/

use dataset1960, clear
egen geoidhhwt=sum(hhwt), by(geoid statefip)
merge m:1 year using choiceadj
keep if _merge==3
drop _merge
replace choiceg=choiceg*adjgas
replace choicee=choicee*adjelec
replace choiceo=choiceo*adjoil
replace choicep=choicep*adjprop
gen indeleccons=hhwt*(predcomeleccons+predelecwatercons*choicewatere+predelecspacecons*choicee)/geoidhhwt
gen indgascons=hhwt*(prednatgaswatercons*choicewaterg+prednatgasspacecons*choiceg)/geoidhhwt
gen indoilcons=hhwt*(predoilwatercons*choicewatero+predoilspacecons*choiceo)/geoidhhwt
gen indpropanecons=hhwt*(predpropanewatercons*choicewaterp+predpropanespacecons*choicep)/geoidhhwt
collapse (sum) indeleccons indgascons indoilcons indpropanecons, by(year statefip geoid)
merge 1:m geoid statefip year using geoid2cty
keep if _merge==3
drop _merge
save countyresults, replace


forvalues y=1961/2019{
use dataset`y', clear
merge m:1 year using choiceadj
keep if _merge==3
drop _merge
replace choiceg=choiceg*adjgas
replace choicee=choicee*adjelec
replace choiceo=choiceo*adjoil
replace choicep=choicep*adjprop
egen geoidhhwt=sum(hhwt), by(geoid statefip)
gen indeleccons=hhwt*(predcomeleccons+predelecwatercons*choicewatere+predelecspacecons*choicee)/geoidhhwt
gen indgascons=hhwt*(prednatgaswatercons*choicewaterg+prednatgasspacecons*choiceg)/geoidhhwt
gen indoilcons=hhwt*(predoilwatercons*choicewatero+predoilspacecons*choiceo)/geoidhhwt
gen indpropanecons=hhwt*(predpropanewatercons*choicewaterp+predpropanespacecons*choicep)/geoidhhwt
collapse (sum) indeleccons indgascons indoilcons indpropanecons, by(year statefip geoid)
merge 1:m geoid statefip year using geoid2cty
keep if _merge==3
drop _merge
append using countyresults
save countyresults, replace
}

/* Step 6: Adjust for households */
use countyresults, clear
replace ctyfip=86 if ctyfip==25 & statefip==12
merge 1:1 statefip ctyfip year using hhnum
keep if _merge==3
drop _merge geoid hhgrowth

local varlist elec gas oil propane
foreach i of local varlist{
gen cty`i'cons=hhnum*ind`i'cons
}
save countyresults, replace

use countyresults, clear
collapse (sum) stateeleccons=ctyeleccons stateoilcons=ctyoilcons stategascons=ctygascons statepropanecons=ctypropanecons, by(statefip year)
format stateeleccons stateoilcons stategascons statepropanecons %12.0f
local varlist gas oil propane elec
foreach i of local varlist{
replace state`i'cons=state`i'cons/1000000
label var state`i'cons "estimated trillions of Btus of `i' consumed"
}

merge 1:1 year statefip using actualenergycons
drop _merge


local varlist gas oil elec propane
foreach i of local varlist{
gen state`i'diff=actstate`i'cons-state`i'cons
}

save stateresults, replace

use stateresults, clear
merge m:1 statefip using stateid
keep if _merge==3
drop _merge state region stateid
collapse (sum) regeleccons=stateeleccons regoilcons=stateoilcons reggascons=stategascons regpropanecons=statepropanecons actreggascons=actstategascons actregoilcons=actstateoilcons actregpropanecons=actstatepropanecons actregeleccons=actstateeleccons, by(year regionid)

merge m:1 year regionid using regadjseds
drop _merge

local varlist gas oil propane elec
foreach i of local varlist{
replace actreg`i'cons=actreg`i'cons/regadj`i'
gen reg`i'diff=actreg`i'cons-reg`i'cons
}
gen regtotaldiff=reggasdiff+regoildiff+regelecdiff+regpropanediff
drop regadjgas regadjelec regadjoil regadjpropane
save regresults, replace

use regresults, clear
collapse (sum) useleccons=regeleccons usoilcons=regoilcons usgascons=reggascons uspropanecons=regpropanecons actuseleccons=actregeleccons actusgascons=actreggascons actusoilcons=actregoilcons actuspropanecons=actregpropanecons, by(year)
format useleccons usoilcons usgascons uspropanecons actuseleccons actusgascons actusoilcons actuspropanecons %12.0f

local varlist gas oil propane elec
foreach i of local varlist{
gen us`i'diff=actus`i'cons-us`i'cons
}
gen ustotaldiff=usgasdiff+usoildiff+uselecdiff+uspropanediff
save usresults, replace






/*Translate to county-level results*/






/*Step 1: Adjust from region to state*/
use actualenergycons, clear
merge m:1 statefip using stateid
keep if _merge==3
drop _merge state region stateid

local varlist gas elec oil propane
foreach i of local varlist{
egen actreg`i'cons=sum(actstate`i'cons), by(year regionid)
gen state`i'share=actstate`i'cons/actreg`i'cons
}

keep year statefip regionid stategasshare stateelecshare stateoilshare statepropaneshare

merge m:1 year regionid using regresults
drop _merge

local varlist gas elec oil propane
foreach i of local varlist{
gen actstate`i'cons=state`i'share*actreg`i'cons
gen state`i'cons=state`i'share*reg`i'cons
}

keep year statefip regionid actstategascons actstateeleccons actstateoilcons actstatepropanecons stategascons stateeleccons stateoilcons statepropanecons

save resultsall, replace


/*Step 2: Create county totals by population*/

use hhnum, clear
drop if year<1960
egen statehh=sum(hhnum), by(statefip year)
gen ctyshare=hhnum/statehh
keep ctyshare year statefip ctyfip
merge m:1 year statefip using resultsall
drop _merge

local varlist gas elec oil propane
foreach i of local varlist{
gen actcty`i'cons=ctyshare*actstate`i'cons
gen cty`i'cons=ctyshare*state`i'cons
}

keep year statefip ctyfip regionid actctygascons actctyoilcons actctyeleccons actctypropanecons ctygascons ctyoilcons ctyeleccons ctypropanecons
save resultsall, replace

/*Step 3: Adjust county totals for consumption differences*/

use countyresults, clear
merge 1:1 year statefip ctyfip using hhnum
keep if _merge==3
drop _merge
egen statehh=sum(hhnum), by(statefip year)

local varlist gas elec oil propane
foreach i of local varlist{
gen ctyshare`i'=hhnum*ind`i'cons/statehh
egen avgstate`i'cons=sum(ctyshare`i'), by(year statefip)
gen ctyadj`i'=ind`i'cons/avgstate`i'cons
}

keep statefip year ctyfip ctyadjgas ctyadjoil ctyadjelec ctyadjpropane
merge 1:1 year statefip ctyfip using resultsall
drop _merge

local varlist gas elec oil propane
foreach i of local varlist{
replace actcty`i'cons=ctyadj`i'*actcty`i'cons
replace cty`i'cons=ctyadj`i'*cty`i'cons
gen cty`i'diff=actcty`i'cons-cty`i'cons
}

save resultsall, replace



/*Create regional and national results*/

use resultsall, clear
replace ctyfip=102 if ctyfip==113 & statefip==46
merge 1:1 statefip ctyfip year using envdamage
keep if _merge==3
gen ctygasdamdiff=ngdam*ctygasdiff
gen ctyoildamdiff=hodam*ctyoildiff
gen ctypropdamdiff=propanedam*ctypropanediff
gen ctyelecdamdiff=elecdam*ctyelecdiff
replace ctyelecdamdiff=0 if ctyelecdamdiff==.

collapse (sum) regiongasdamdiff=ctygasdamdiff regionoildamdiff=ctyoildamdiff regionpropdamdiff=ctypropdamdiff regionelecdamdiff=ctyelecdamdiff regionactgascons=actctygascons regionactoilcons=actctyoilcons regionacteleccons=actctyeleccons regionactpropcons=actctypropanecons regiongascons=ctygascons regionoilcons=ctyoilcons regioneleccons=ctyeleccons regionpropcons=ctypropanecons, by(region year)
save regionresults, replace

collapse (sum) usgasdamdiff=regiongasdamdiff usoildamdiff=regionoildamdiff uspropdamdiff=regionpropdamdiff uselecdamdiff=regionelecdamdiff usactgascons=regionactgascons usactoilcons=regionactoilcons usactpropcons=regionactpropcons usacteleccons=regionacteleccons usgascons=regiongascons usoilcons=regionoilcons uspropcons=regionpropcons useleccons=regioneleccons, by(year)
save usresults, replace






/*Tables and Figures*/






/*Table 2*/

cd "C:\Users\amhill3\Desktop\Projects\old or done\Published\NatGas Deregulation Env\Data\estimation"
use fuel1, clear
keep fuelheat fuelwatr regionc hhwt year
gen gash=0
replace gash=hhwt if fuelheat==2
gen proph=0
replace proph=hhwt if fuelheat==3
gen elech=0
replace elech=hhwt if fuelheat==4
gen oilh=0
replace oilh=hhwt if fuelheat==5

gen gasw=0
replace gasw=hhwt if fuelwatr==2
gen propw=0
replace propw=hhwt if fuelwatr==3
gen elecw=0
replace elecw=hhwt if fuelwatr==4
gen oilw=0
replace oilw=hhwt if fuelwatr==5

egen usgash=sum(gash), by(year)
egen usoilh=sum(oilh), by(year)
egen usproph=sum(proph), by(year)
egen uselech=sum(elech), by(year)
egen ushhwt=sum(hhwt), by(year)
egen reggash=sum(gash), by(year regionc)
egen regoilh=sum(oilh), by(year regionc)
egen regproph=sum(proph), by(year regionc)
egen regelech=sum(elech), by(year regionc)
egen reghhwt=sum(hhwt), by(year regionc)

egen usgasw=sum(gasw), by(year)
egen usoilw=sum(oilw), by(year)
egen uspropw=sum(propw), by(year)
egen uselecw=sum(elecw), by(year)
egen reggasw=sum(gasw), by(year regionc)
egen regoilw=sum(oilw), by(year regionc)
egen regpropw=sum(propw), by(year regionc)
egen regelecw=sum(elecw), by(year regionc)

gen pusgash=usgash/ushhwt
gen pusoilh=usoilh/ushhwt
gen pusproph=usproph/ushhwt
gen puselech=uselech/ushhwt

gen preggash=reggash/reghhwt
gen pregoilh=regoilh/reghhwt
gen pregproph=regproph/reghhwt
gen pregelech=regelech/reghhwt

gen pusgasw=usgasw/ushhwt
gen pusoilw=usoilw/ushhwt
gen puspropw=uspropw/ushhwt
gen puselecw=usgasw/ushhwt

gen preggasw=reggasw/reghhwt
gen pregoilw=regoilw/reghhwt
gen pregpropw=regpropw/reghhwt
gen pregelecw=regelecw/reghhwt

collapse pusgash pusoilh pusproph puselech preggash pregoilh pregproph pregelech pusgasw pusoilw puspropw puselecw preggasw pregoilw pregpropw pregelecw, by(regionc year)

/*Get 1990 water heating from 1990 RECS*/
use RECSadj, clear
keep if year==1990
keep fuelh2o regionc nweight year

gen gasw=0
replace gasw=nweight if fuelh2o==1
gen propw=0
replace propw=nweight if fuelh2o==2
gen oilw=0
replace oilw=nweight if fuelh2o==3
gen elecw=0
replace elecw=nweight if fuelh2o==5

collapse (sum) gasw propw oilw elecw nweight, by(regionc year)
reshape wide gasw propw oilw elecw nweight, i(year) j(regionc)

local varlist gasw propw oilw elecw nweight
foreach i of local varlist{
gen `i'0=`i'1+`i'2+`i'3+`i'4
}
reshape long gasw propw oilw elecw nweight, i(year) j(regionc)

replace gasw=gasw/nweight
replace oilw=oilw/nweight
replace propw=propw/nweight
replace elecw=elecw/nweight
drop nweight


/*Figure 1*/


use hddgraphic, clear
egen regpop=sum(pop), by(region year)
gen regshare=hdd*pop/regpop
egen reghdd=sum(regshare), by(region year)
collapse reghdd, by(region year)
reshape wide reghdd, i(year) j(region)


/*Table 3*/


use pricegraphic, clear
egen regpop=sum(pop), by(region year)
gen ygroup=1
replace ygroup=2 if year>1973
replace ygroup=3 if year>1986
replace ygroup=4 if year>2004
local varlist natgas propane elec oil
foreach i of local varlist{
gen regshare`i'=p`i'*pop/regpop
egen regp`i'=sum(regshare`i'), by(year region)
}

collapse regpnatgas regppropane regpelec regpoil, by(ygroup region)
keep regpnatgas regppropane regpelec regpoil ygroup region
reshape long regp, i(ygroup region) j(fuel) string
reshape wide regp, i(fuel ygroup) j(region)


/*Figure 2*/


use ctyemissions, clear
replace ctyfip=86 if ctyfip==25 & statefip==12
merge m:1 year statefip ctyfip using damfac
keep if _merge==3
gen pm25damtot=pm25dam*pm25tons
gen so2damtot=so2dam*so2tons
gen co2damtot=co2dam*co2tons
gen noxdamtot=noxdam*noxtons
collapse (sum) pm25damtot so2damtot co2damtot noxdamtot, by(year nerc)
merge 1:1 year nerc using nercmwh70-11
drop _merge
gen dampm25=pm25damtot/mwh
gen damco2=co2damtot/mwh
gen damso2=so2damtot/mwh
gen damnox=noxdamtot/mwh
gen elecdam=dampm25+damco2+damso2+damnox
collapse elecdam, by(year nerc)
reshape wide elecdam, i(year) j(nerc) string
set obs 43
replace year=2001 in 38
replace year=2002 in 39
replace year=2003 in 40
replace year=2006 in 41
replace year=2008 in 42
replace year=2011 in 43
sort year
local varlist ECAR ERCOT FRCC MAAC MAIN MAPP NPCC SERC SPP WECC
forvalues y=2001/2003{
foreach i of local varlist{
replace elecdam`i'=elecdam`i'[_n-1] if year==`y'
}
}

local varlist ERCOT FRCC MRO NPCC RFC SERC SPP WECC
foreach i of local varlist{
replace elecdam`i'=(elecdam`i'[_n-1]+elecdam`i'[_n+1])/2 if year==2006
replace elecdam`i'=(elecdam`i'[_n-1]+elecdam`i'[_n+1])/2 if year==2008
replace elecdam`i'=(elecdam`i'[_n-1]+elecdam`i'[_n+1])/2 if year==2011
}
reshape long elecdam, i(year) j(nerc) string
reshape wide elecdam, i(nerc) j(year)



/*Create national results without env damage*/

use resultsall, clear
collapse (sum) actctygascons ctygascons actctyeleccons ctyeleccons actctyoilcons ctyoilcons actctypropanecons ctypropanecons, by(year)
gen diffgas=actctygascons-ctygascons
gen diffoil=actctyoilcons-ctyoilcons
gen diffelec=actctyeleccons-ctyeleccons
gen diffpropane=actctypropanecons-ctypropanecons
keep diffgas diffoil diffelec diffpropane year
save Fig3, replace

/*Gas and Oil focus*/
use resultsall, clear
collapse (sum) actctygascons ctygascons actctyoilcons ctyoilcons, by(year)
save Fig4, replace

/*Create regional results without env damage*/
use resultsall, clear
collapse (sum) actctygascons ctygascons actctyeleccons ctyeleccons actctyoilcons ctyoilcons actctypropanecons ctypropanecons, by(regionid year)
gen diffgas=actctygascons-ctygascons
gen diffoil=actctyoilcons-ctyoilcons
gen diffelec=actctyeleccons-ctyeleccons
gen diffpropane=actctypropanecons-ctypropanecons
keep regionid year diffgas diffoil diffelec diffpropane
save Fig5, replace

/*Gas shortage by region*/
use usresults, clear
keep usactgascons usgascons year
gen gasshortage=-(usactgascons-usgascons)/usgascons
gen period=1
replace period=2 if year>1969
replace period=3 if year>1979
replace period=4 if year>1989
replace period=5 if year>1999
replace period=6 if year>2009
collapse gasshortage, by(period)
gen regionid=5
save Tab6, replace

use resultsall, clear
collapse (sum) actctygascons ctygascons, by(regionid year)
gen gasshortage=-(actctygascons-ctygascons)/ctygascons
gen period=1
replace period=2 if year>1969
replace period=3 if year>1979
replace period=4 if year>1989
replace period=5 if year>1999
replace period=6 if year>2009
collapse gasshortage, by(period regionid)
append using Tab6
save Tab6, replace

/*Create placebo test*/
use stateresults, clear
keep if statefip==22 | statefip==40 | statefip==35 | statefip==48
collapse (sum) actstategascons stategascons, by(year)
save Fig6, replace

/*Create emissions difference*/
use resultsall, clear

merge 1:1 year statefip ctyfip using elecemissions
drop if _merge==2
replace elecpm25=0 if _merge==1
replace elecco2=0 if _merge==1
replace elecso2=0 if _merge==1
replace elecnox=0 if _merge==1
keep elecpm25 elecco2 elecso2 elecnox statefip year ctyfip regionid ctyoildiff ctyelecdiff ctypropanediff ctygasdiff
merge m:1 year using emissionscontentbyfuel
keep if _merge==3
drop _merge
rename ctyoildiff ctyhodiff
rename ctygasdiff ctyngdiff

local varlist ho ng elec propane
local varlist2 co2 pm25 so2 nox
foreach i of local varlist{
foreach j of local varlist2{
gen `j'`i'=cty`i'diff*`i'`j'
}
}
foreach j of local varlist2{
gen `j'tot=1000*(`j'ng+`j'ho+`j'elec+`j'propane)
}

collapse (sum) co2tot so2tot noxtot pm25tot, by(year regionid)
save emissionsimpact, replace
collapse (sum) co2tot so2tot noxtot pm25tot, by(year)
gen regionid=5
append using emissionsimpact
gen period=1
replace period=2 if year>1984
collapse co2tot so2tot noxtot pm25tot, by(regionid period)
save Tab7, replace

/*National Environmental Damage*/
use usresults, clear
gen totusdam=(usgasdamdiff+usoildamdiff+uspropdamdiff+uselecdamdiff)/1000
keep totusdam year
save Fig7, replace

/*Regional Environmental Damage*/
use regionresults, clear
gen totregiondam=(regiongasdamdiff+regionoildamdiff+regionpropdamdiff+regionelecdamdiff)/1000
keep totregiondam year regionid
gen period=1
replace period=2 if year>1969
replace period=3 if year>1979
replace period=4 if year>1989
replace period=5 if year>1999
replace period=6 if year>2009
collapse totregiondam, by(regionid period)
reshape wide totregiondam, i(period) j(regionid)
gen totregiondam5=totregiondam1+totregiondam2+totregiondam3+totregiondam4
save Tab8, replace

/*Create emissions damage by fuel source*/
use usresults, clear
gen totdam=(usgasdamdiff+usoildamdiff+uspropdamdiff+uselecdamdiff)/1000
gen gasdamdiff=usgasdamdiff/1000
gen oildamdiff=usoildamdiff/1000
gen propdamdiff=uspropdamdiff/1000
gen elecdamdiff=uselecdamdiff/1000
keep totdam gasdamdiff oildamdiff propdamdiff elecdamdiff
collapse totdam gasdamdiff oildamdiff propdamdiff elecdamdiff
gen regionid=5
save Tab9, replace

use regionresults, clear
gen totdam=(regiongasdamdiff+regionoildamdiff+regionpropdamdiff+regionelecdamdiff)/1000
gen gasdamdiff=regiongasdamdiff/1000
gen oildamdiff=regionoildamdiff/1000
gen propdamdiff=regionpropdamdiff/1000
gen elecdamdiff=regionelecdamdiff/1000
keep totdam gasdamdiff oildamdiff propdamdiff elecdamdiff regionid
collapse totdam gasdamdiff oildamdiff propdamdiff elecdamdiff, by(regionid)
append using Tab9
save Tab9, replace

/*Create county results for map*/

use resultsall, clear
replace ctyfip=102 if ctyfip==113 & statefip==46

merge 1:1 statefip ctyfip year using envdamage
keep if _merge==3
gen ctygasdamdiff=ngdam*ctygasdiff
gen ctyoildamdiff=hodam*ctyoildiff
gen ctypropdamdiff=propanedam*ctypropanediff
gen ctyelecdamdiff=elecdam*ctyelecdiff
replace ctyelecdamdiff=0 if ctyelecdamdiff==.
gen totdiff=ctygasdamdiff+ctyoildamdiff+ctypropdamdiff+ctyelecdamdiff
keep statefip year ctyfip totdiff
collapse (sum) totdiff, by(statefip ctyfip)
save Fig8, replace

/*Top 10 counties*/

use resultsall, clear
replace ctyfip=102 if ctyfip==113 & statefip==46

merge 1:1 statefip ctyfip year using envdamage
keep if _merge==3
gen ctygasdamdiff=ngdam*ctygasdiff
gen ctyoildamdiff=hodam*ctyoildiff
gen ctypropdamdiff=propanedam*ctypropanediff
gen ctyelecdamdiff=elecdam*ctyelecdiff
replace ctyelecdamdiff=0 if ctyelecdamdiff==.
gen totdiff=ctygasdamdiff+ctyoildamdiff+ctypropdamdiff+ctyelecdamdiff
keep statefip year ctyfip totdiff year
gen id=0
replace id=1 if statefip==36 & ctyfip==81
replace id=1 if statefip==34 & ctyfip==3
replace id=1 if statefip==25 & ctyfip==17
replace id=1 if statefip==36 & ctyfip==47
replace id=1 if statefip==34 & ctyfip==13
replace id=1 if statefip==36 & ctyfip==59
replace id=1 if statefip==42 & ctyfip==3
replace id=1 if statefip==17 & ctyfip==31
replace id=1 if statefip==25 & ctyfip==27
replace id=1 if statefip==26 & ctyfip==163
keep if id==1
gen period=1
replace period=2 if year>1984
collapse totdiff, by(statefip ctyfip period)
reshape wide totdiff, i(statefip ctyfip) j(period)
sort statefip ctyfip
save Tab10, replace